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ABSTRACT 

In the magnetar model, the quiescent non-thermal soft X-ray emission from Anomalous X-ray Pulsars and 
Soft-Gamma Repeaters is thought to arise from resonant comptonization of thermal photons by charges mov- 
ing in a twisted magnetosphere. Robust inference of physical quantities from observations is difficult, because 
the process depends strongly on geometry and current understanding of the magnetosphere is not very deep. 
The polarization of soft X-ray photons is an independent source of information, and its magnetospheric imprint 
remains only partially explored. In this paper we calculate how resonant cyclotron scattering would modify the 
observed polarization signal relative to the surface emission, using a multidimensional Monte Carlo radiative 
transfer code that accounts for the gradual coupling of polarization eigenmodes as photons leave the magne- 
tosphere. We employ a globally-twisted, self-similar, force-free magnetosphere with a power-law momentum 
distribution, assume a blackbody spectrum for the seed photons, account for general relativistic light deflection 
close to the star, and assume that vacuum polarization dominates the dielectric properties of the magnetosphere. 
The latter is a good approximation if the pair multiplicity is not much larger than unity. Phase-averaged po- 
larimetry is able to provide a clear signature of the magnetospheric reprocessing of thermal photons and to 
constrain mechanisms generating the thermal emission. Phase-resolved polarimetry, in addition, can charac- 
terize the spatial extent and magnitude of the magnetospheric twist angle at ^ 100 stellar radii, and discern 
between uni- or bidirectional particle energy distributions, almost independently of every other parameter in 
the system. We discuss prospects for detectability with GEMS. 
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1. INTRODUCTION 

The rich phenomenology associated with Soft Gamma Re- 
peaters (SGRs; e.g. Mereghetti et al. 2007) and Anomalous 
X-ray Pulsars (AXPs; e.g. Kaspi 2007) has so far been best 
accounted for by the magnetar model, which hypothesizes the 
decay of a strong magnetic field in a neutron star as the dom- 
inant energy source (Duncan & Thompson 1992; Thompson 
& Duncan 1996). 

In the soft X-ray band (1-10 keV), the quiescent spectra of 
magnetar candidates are, to first approximation, very similar: 
a thermal component with temperature kT ^ QA keV, and a 
non-thermal tail of varying hardness at higher energies (e.g.. 
Woods & Thompson 2006). On an individual basis, though, 
quiescent AXPs and SGRs can show variability in their spec- 
tra and pulse profiles over timescales ranging from days to 
years, in many cases around outburst periods or glitches (e.g., 
Kaspi 2007; Woods et al. 2007). 

In the magnetar model, the non-thermal component of the 
spectrum arises from the decay of the toroidal part of the mag- 
netospheric field (the "twist"), thought to be acquired during 
outbursts (Thompson et al. 2002). The twist induces a current 
along the field lines, increasing the density of charges by or- 
ders of magnitude above the Goldreich & Julian (1969) value, 
resulting in an optical depth to resonant cyclotron scattering of 
order unity. Comptonization of thermal photons from the stel- 
lar surface then generates a non-thermal tail in the spectrum, 
and pulsed emission from the combination of stellar rotation 
and the geometric dependence of the process (Thompson et al. 
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2002). 

Numerical calculations of the resulting emission have been 
carried out with both one-dimensional semianalytic methods 
(Lyutikov & Gavriil 2006) and multidimensional Monte Carlo 
algorithms (Fernandez & Thompson 2007; Nobili et al. 2008). 
Output from the one-dimensional computation has been used 
to model phase averaged spectra of AXPs and SGRs (Giiver 
et al. 2007, 2008; Rea et al. 2007a, 2008), resulting in good 
fits while improving upon the ad-hoc blackbody plus power- 
law prescription. Even better fits are obtained when using 
the results of the Nobili et al. (2008) calculation, as no extra 
power-law is required to fit SGR spectra (Zane et al. 2009). 

Still, phase averaged observables convey only part of the 
information, raising the question of geometric degeneracies. 
Recently, Albano et al. (2010) have performed simultaneous 
light curve and spectral fitting of the transient AXPs XTE 
J18 10-197 and CXOU J 1647 10.2-4552 16, making use of the 
Nobili et al. (2008) output, constraining for the first time 
the geometric orientation of these sources and hence provid- 
ing a more robust spectral fit. Nevertheless, both the Fer- 
nandez & Thompson (2007) and Nobili et al. (2008) studies 
make reasonable but heuristic assumptions about the magne- 
tospheric physics, caused largely by the difficulty of modeling 
the closed field line circuit with radiation feedback from first 
principles. 

Barring an improved understanding of the physics, fur- 
ther progress is possible by making use of information en- 
coded in the polarization signal. In recent years, photoelec- 
tric polarimeters have reached the sensitivity level required 
for detecting polarized X-rays from a variety of astrophysical 
sources (Costa et al. 2001; Black et al. 2007), allowing exten- 
sion of the discovery space beyond the lone data point for the 
Crab nebula (Novick et al. 1972). In particular, the Gravity 
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and Extreme Magnetism (GEMS) mission will reach a Mini- 
mum Detectable Polarization (MDP) of ^ 1 % for a 1 mCrab 
source in the 2- 10 keV range with a 10* s integration (Swank 
et al. 2010). About one half of the current magnetar candi- 
dates are brighter than 1 mCrab'', making detection of their 
polarization a concrete possibility. 

Polarized X-ray transfer calculations have been carried 
out for isolated neutron stars with dipole-type magnetic 
fields, in some cases including realistic atmospheres but with 
no magnetospheric processes other than propagation in a 
Schwarzschild spacetime (Heyl & Shaviv 2000; Pavlov & 
Zavlin 2000; Heyl et al. 2003; van Adelsberg & Lai 2006; 
Wang & Lai 2007; van Adelsberg & Perna 2009; Wang & Lai 
2009). Both the Fernandez & Thompson (2007) and NobiU 
et al. (2008) Monte Carlo studies assume adiabatic propaga- 
tion, with photons remaining in a fixed polarization eigen- 
mode in between scatterings, until they escape the magne- 
tosphere. Nobili et al. (2008) estimated the observed polar- 
ization fractions through the relative difference between the 
emerging normal mode intensities. The polarization measured 
by a distant observer, however, is set by a gradual coupling of 
the propagation eigenmodes as photons leave the magneto- 
sphere, with the generation of a circular component in some 
cases (Heyl & Shaviv 2000). Thus, evolution of all polariza- 
tion degrees of freedom is required for making observational 
predictions. 

In this paper we extend the multidimensional Monte Carlo 
radiative transfer code of Fernandez & Thompson (2007) to 
include polarization propagation and general relativistic light 
deflection close to the neutron star. Our focus is to under- 
stand the basic features imparted to the polarization signal by 
the scattering process and the geometry of the system. To 
this end, we perform radiative transfer simulations using the 
simplest possible prescriptions for magnetic field geometry, 
particle energy distribution, and seed photons, when the di- 
electric properties of the medium are dominated by vacuum 
polarization. Motivated by GEMS, we explore the parameter 
dependencies of polarization fraction and polarization angle 

The structure of this paper is the foUowing. Section 2 lays 
out the physics included in the model, while Section 3 dis- 
cusses our numerical implementation. Section 4 contains our 
results, focusing on the dependence of observables on sys- 
tem parameters. In Section 5 we discuss our findings in light 
of the expected instrumental parameters of GEMS, and esti- 
mate integration times for a few bright magnetar candidates. 
A summary and discussion of our results is provided in Sec- 
tion 6; readers not interested in technical details can start here, 
where explicit reference is made to the key figures and equa- 
tions. Appendix A contains details about geometric transfor- 
mations employed in the code. 

2. PHYSICAL MODEL 

The elements influencing the polarization signal measured 
at infinity are: (i) the physics of the magnetosphere, includ- 
ing the transfer of polarized radiation, (ii) the properties of the 
seed photons, and (iii) the type of spacetime. We describe our 
treatment of each of these items in turn. A detailed description 
of the physics entering the twisted magnetosphere model can 
be found in Thompson et al. (2002) and Fernandez & Thomp- 
son (2007). Here we provide a minimal account, focused on 
the part relevant to polarization observables. 

* A catalog of current magnetar candidates is maintained at 
http://www.physics.mcgill.ca/~pulsar/magnetar/main.html 



2.1. Twisted Magnetosphere 

The basic hypothesis behind the magnetar model is the ex- 
istence of a strong toroidal field inside the star, thought to 
have been generated via dynamo action in a rapidly rotating 
protoneutron star (Thompson & Duncan 1995; Thompson & 
Duncan 2001). Seismic instabihties in the crust induced by 
strong magnetic stresses mediate the transfer of some of this 
internal toroidal energy to the magnetosphere, imparting a fi- 
nite twist to the otherwise dipolar external field (Thompson & 
Duncan 200 1 ; Thompson et al. 2002; Beloborodov & Thomp- 
son 2007). 

To describe such a magnetosphere, one would ideaUy solve 
a closed field line circuit including pair creation processes 
and feedback from radiation at cyclotron frequencies (e.g., 
Thompson et al. 2002; Beloborodov & Thompson 2007; 
Thompson 2008a,b). The global problem has not yet been 
solved self-consistently. Numerical computations of the X- 
ray emission rely instead on prescriptions for the field geom- 
etry and particle energy distribution that can capture the basic 
aspects of the model (Fernandez & Thompson 2007; Nobili 
et al. 2008; Pavan et al. 2009). As our goal here is to explore 
the dependence of the polarization signal on different param- 
eters rather than to fit individual sources, we adopt the same 
parametric approach in our calculations. 

2.1.1. Magnetic Field Geometry and Particle Energy 
Distribution 

We adopt the force-free, self-similar, twisted-dipole solu- 
tion of Thompson et al. (2002), 

1 f R \ 
B = -Bpole(^j F(cos^), (1) 

with r the radial distance from the center of the star, 9 the po- 
lar angle relative to the magnetic axis, Bpoie the strength of the 
field at the poles, R^s the stellar radius, and p a constant. The 
components of the function F are Fr = -/', F$ = {p/ sm6)f, 
and = [C/{p(p+ l)}]f^''Fe, with /(cos6i) the solution to 
the second order differential equation 



sm^ef" + Cf^^/P + p(p+l)f=0, 



(2) 



subject to the boundary conditions /'(O) = 0, /'(I) = -2, and 
/(I) = 0. This leaves an additional free parameter, either C or 
p. Alternatively, one can label solutions by the net twist angle 
of the field Unes that are anchored close to the magnetic poles. 
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When going from A0ai_s = to tt, the solutions interpolate 
smoothly between a dipole (p = I) and a split monopole (p = 
0). For illustration, Figure 1 shows a few field Unes anchored 
at 15° from the magnetic axis for a twist A(^n-s = 1- 

The current density induced along the twisted hnes is 
(Thompson et al. 2002) 



J = y^Z,en,/3,Bc = 



(p+l)cB4 
Airr Bg 



B, 



(4) 



where Z,e is the electric charge of species i, n,- the number 
density, B = B/B, and 



% = J fi(r,p)^dp 



(5) 
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Fig. 1 . — Diagram showing tlie spatial an'angement of tlie problem (not to 
scale). The neutron star (grey circle) emits thermal photons (green trajectory) 
which scatter electrons moving along field lines (black solid lines). Magne- 
tospheric charges are provided by a current J due to the net magnetospheric 
twist (§2.1.1). Photons can scatter only inside r^sc (eq. [26]). At radii much 
smaller than Tpi (eq. [27]), the polarization propagation is adiabatic (dashed 
green line, §2.1.3). When approaching this radius, polarization eigenmodes 
start to couple (solid green) and eventually freeze beyond r^i (dot-dashed 
green). For clarity, only a few field lines anchored at 15° from the magnetic 
axis are shown, for a twist angle A(}!)im_s = 1 (eq. [3]). Photon trajectories are 
curved close to the star due to light-bending (§2.3). The vectors M and Q 
denote the magnetic and rotation axes, respectively. 



the mean velocity along the field lines in units of c, with 
/(r , p) the one-dimensional phase-space distribution function, 
normalized such that J f(r,p)dp = 1. The number density of 
charges n, is obtained from 



|Z,|e«,_ (p+l)B4 



B 



A-KrPi Be 



(6) 



where e, is the fraction of the current carried by species i. 

To complete the description of the magnetosphere one 
needs to specify the particle energy distribution of the charge 
carriers. We employ the simplest possible prescription, a 
power law in momentum, independent of position. 



/(7/?)cx(7/?r", 



(7) 



where 7 = (1 - /3^)"'/^ is the Lorentz factor. This distribution is 
characterized by three parameters: the exponent a, the mini- 
mum velocity /S^in, and the maximum Lorentz factor 7max- In 
all of our runs we employ a = -2 and /3min = 0.2, but allow 
7max to vary (see, e.g., Fernandez & Thompson 2007 for the 
resulting spectra and their comparison with other distribution 
functions). To account for both an electron-ion and electron- 
positron plasma in the magnetosphere, we allow equation (7) 
to be uni- or bi-directional, that is, extending over only posi- 
tive (from north to south magnetic pole. Fig. [1]) or both pos- 
itive and negative momenta^. 

^ In the case of an electron-ion plasma, the relevant charge carrier for the 
interaction with keV X-rays far from the star are only the electrons. 



In the simple case considered here {tie ~ -//[ec]), the di- 
electric and inverse permeability tensors are dominated by 
vacuum polarization (§2.1.3). As a first step in addressing 
this problem, and to ease comparison with previous work, 
we ignore the possibility of a large pair multiplicity, in which 
case lie 3> J /{ec). Such an effect could result from pair cas- 
cades originated by the energy deposition of current-driven in- 
stabilities and the rapid conversion of resonantly upscattered 
X-rays into pairs under a strong magnetic field (Thompson 
2008a). This could significantly modify the dielectric proper- 
ties at distances from the star where the adiabatic approxima- 
tion breaks down. Given that our modeling of the magneto- 
spheric currents is rather heuristic, the inclusion of plasma po- 
larization could lead to artificial features sensitive to the tem- 
perature and density of the magnetospheric particles. Hence, 
we opt for a self-consistent but potentially incomplete calcu- 
lation, and only include vacuum polarization. Consequently, 
we set e, = e = 1 in equation (6) for an ion-electron plasma, or 
e,- =1/2 when electrons and positrons are present. 

2.1.2. Resonant Cyclotron Scattering 

For magnetar field strengths, the lifetime of 
an excited Landau level is very short, AfL = 
(3/4)(rV[a™m,c2])(BQED/B)2 3 x lQ-^%\ s, with a,^ the 
fine sti-ucture constant, Bqed = meC^/(He) ~ 4.4 x lO'-' G the 
field strength where the cyclotron energy equals the electron 
rest mass energy, and Bn =B/(10"G). Hence, electrons 
or positrons can be assumed to be in the ground state, and 
absorption plus re-emission of photons can be treated as a 
single scattering process (Zheleznyakov 1996). Excitation 
occurs when the photon frequency, in the rest frame of the 
charge, equals the cyclotron frequency. In the stellar frame, 
this condition corresponds to 



7(1 -/3m) 

where cj is the photon angular frequency, 
p = k-B 



(8) 



(9) 



is the cosine of the angle between photon (k) and magnetic 
field directions in the stellar frame, and 



UJr= ■ 



eB 
nigC 
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(10) 



the electron cyclotron frequency. 

The total cross section in the stellar frame is given by (e.g., 
Meszaros 1992) 



is 



(11) 



where \e±\ is the rest frame overlap of the photon polariza- 
tion state with a left (-) or right (-I-) circularly polarized wave, 
depending on the sign of the charge. The factor (l-/3/i) 
comes from the relativistic boost to the cross section (e.g., 
Rybicki & Lightman 2004). The delta function is an ap- 
proximation to the small fractional width of the resonance for 
a single particle, (lAti^cUeT^ = (2a^j3)(B/BQED) = IQ-^Bu- 
The enhancement relative to the usual non-resonant Thomson 
cross section ctj- is, at resonance and in the charge rest frame, 
a,^,/cTT = |e±p(2AfLWc)^ - 10'" |e±pB?i (Meszaros 1992). 

Calculation of the optical depth requires averaging the cross 
section over the phase space distribution of the scattering 
charges. Using the fact that for fixed photon frequency and 
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at a given position, equation (8) has two roots /3^, the differ- 
ential optical depth along a small distance A£ <C f is given by 
(Fernandez & Thompson 2007) 



7r(/7+l)w / B 
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Be 
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/(7r)A(7/?)- 



7/3- J 



/(7/3^)A(7/3)+ 
(12) 



where (7/3)^ are the resonant momenta corresponding to the 
resonant velocities (3^, A{jl3)^ is the size of the step in 
resonant-momentum space corresponding to the spatial step 
A£, and dojo/di = k ■ Wujd- Equation (12) differs from Fer- 
nandez & Thompson (2007) in that the average velocity |^| 
comes out of the integral, as it is given by the relation be- 
tween current density and number density (eq. [4]) and does 
not depend on the resonance condition^\ Given that typical 
values of /? are not too far from the mean of the distribution, 
this modification does not introduce a significant difference in 
the results relative to those of Fernandez & Thompson (2007). 

During scattering, a photon can change its polarization 
state. The probabilities depend on the rest frame overlap of 
the outgoing photon polarization with a circularly polarized 
wave, as the differential cross section for scattering from po- 
larization state A to state B satisfies (Meszaros 1992) 



dcr, 



AO. 



res I |2| / |2 



(13) 



where primes denote outgoing photon quantities. If normal 
modes begin to couple, the simple adiabatic expressions for 
the ingoing photon overlap need to be modified, as described 
in §3.2 (eq. [33]). For the outgoing photon, however, we im- 
pose a normal mode state, regardless of the degree of adi- 
abaticity. As vacuum polarization dominates the dielectric 
properties (§2. 1.3), the electric field eigenvectors are given by 
ce oc(k.x B) and eo oc {ce x k). The overlap is the dot prod- 
uct of the normalized eigenvectors with e± = (i ± iy) / \/2 in a 
frame where B is along z, yielding 



= 1/2, K±|2=(^;)V2, 



where 



Mr 



(14) 



(15) 



is the angle between photon direction and magnetic field in 
the rest frame of the charge. The probabiUty of the photon 
coming out in E-mode is then 



P(E) = 



1 



i+(m;)2 



(16) 



with P(0) = 1 —P(E). In practice, imposing a normal mode 
right after scattering does not affect our results, as the polar- 
ization coupling and scattering regions are well separated in 
space (§2.1.3, see also Figure 2). 

Note that for all < 1, P(0) < P(E), and hence E-mode 
photons always dominate the outgoing polarization distribu- 
tion except for exactly parallel propagation. 

* Nobili et al. (2008) used an optical depth with this error corrected. The 
factor ojj^' (duij) / d£) is independent of ojd = c<J, so there is no dependence on 
photon frequency or magnetic field strength. 



2.1.3. Transfer of Polarized Radiation 

We first analyze the dielectric properties of the magne- 
tosphere, then write down the polarization evolution equa- 
tions, and finally characterize the regions over which mag- 
netospheric scattering and non-adiabatic propagation are rel- 
evant, looking for possible overlaps. 

For B <C Bqed, the relative contributions of plasma and 
vacuum polarization to the dielectric and inverse permeabil- 
ity tensors can be quantified by the ratio (e.g., Zheleznyakov 
1996) 



{■uJp/uJ) / 1 A-7 -1 —1(^1 



1+p 



sin^6>A(/)N-s / \ 



where uip = \f^Tx7?ehvJrne is the electron plasma frequency, 
and the current density (eq. [4]) has been approximated by 
J ~ c(47rr)"' sin^^^ A0N-sB (Thompson et al. 2002). Equa- 
tion (17) shows that vacuum polarization dominates out to r ~ 
3000/?NS for mildly relativistic particle distributions (|^| < 1). 
Outside this radius plasma polarization starts to dominate. 

We show in §3.3 that the polarization vector effectively 
freezes at distances < 5007?ns, well inside the vacuum dom- 
inated regime. However, this would change if the pair mul- 
tiplicity is bigger than M.± ^ 10^^, because the radius where 
plasma and vacuum polarization become comparable would 
decrease by a factor (Al±)^/<'+''' 30. 

Due to quantum electrodynamic effects, a strong magnetic 
field modifies the dielectric properties of the vacuum accord- 
ing to (Klein & Nigam 1964; Adler 1971) 



£ = (l+a)I-(-^BB 
/i = (l-l-a)I-l-m5B, 



(18) 
(19) 



where e and ]i are the dielectric and inverse permeability ten- 
sors, respectively, and I is the unit tensor. For B -C Sqed. 
the coefficients are a = -25, q = 15, and m = —A5, with 5 = 
aem(457r)-i(B/BQED)^ ^ 3 X lO-'^B^,. These coefficients are 
all real, so there is no change in the total intensity. The prop- 
agation eigenmodes are linear, and are defined in relation to 
the photon and magnetic field directions (§2.1.2). General 
expressions exist for arbitrary field strengths, which neverthe- 
less maintain the geometric structure of equations (18)-(19) 
[Heyl & Hernquist 1997]. Given that polarization eigenmodes 
begin to decouple at distances r > 10/?ns. the weak field limit 
is good enough for our purposes. 

The polarization of photons can be characterized through 
the complex components of the electric field vector in some 
coordinate system, or through the Stokes parameters (e.g., 
Rybicki & Lightman 2004). In our implementation, we 
choose to evolve the former, as it is simpler to perform ge- 
ometric transformations on them, and convert to Stokes pa- 
rameters only after the polarization is deemed to be frozen 
(3.3). 

We now derive the ordinary differential equations govern- 
ing the evolution of the electric field vector in the geometric 
optics limit, following the notation of Lai & Ho (2003b). The 
wave equation for the electric field of a plane electromagnetic 

wave E = Eo(r)e"''^' is 



CO 



Vx[;U-(VxE)]=^£-E. 



(20) 
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The polarization tensors are assumed to be constant in time 
and varying over distances much larger than the photon wave- 
length. We then choose a coordinate system with the z-axis 
along the direction of photon propagation, assume that quan- 
tities change only along z, and write Eo(z) = exp(jfcoz)A(z), 
with ^0 = w/c. The geometric optics approximation is equiva- 
lent to d^A. <^ koA. Keeping terms linear in k^^dz, we obtain a 
system of ordinary differential equations describing the prop- 
agation of the complex ampUtude A, 

-1 r / 

IJ'XX 
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- H'yx r _ _ 1 
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(22) 



(23) 



Equations (21)-(23) are valid for general conditions, so long 
as the geometric optics and smooth background approxima- 
tion holds. Inclusion of the plasma contribution to fl and e is 
mathematically straightforward, once the number density and 
temperature of charges as a function of position are known 
(see, e.g., Lai & Ho 2003b for an appUcation to neutron star 
atmospheres). 

Even for magnetar field strengths, the non-diagonal com- 
ponents of £ are much smaller that unity, hence from equa- 
tion (23) we can neglect A^. To more clearly visuaUze the 
behavior of equations (21)-(22), one can switch to a ref- 
erence frame where the magnetic field is in the x-z plane, 
B = cos6kBk+sm6kBX, and keep terms Unear in {a,q,m}, find- 
ing 



(24) 



In this case, the x and y components correspond to the O- and 
E-mode eigenvectors, respectively. The characteristic scale 
over which A varies is 

4 = ^, (25) 

where An = sin^ 9kBiq+m) /2 is the difference between the in- 
dices of refraction (Heyl & Shaviv 2000; Lai & Ho 2003b). 
As the photon propagates, non-diagonal components will ap- 
pear in equation (24), because B will leave the x-z plane. 
However, if this occurs over a length-scale much larger than 
Ia, normal modes do not mix, and the propagation is adia- 
batic (Heyl & Shaviv 2000; Lai & Ho 2003b). As the pho- 
ton moves out in the magnetosphere, the QED correction to 
the dielectric tensor falls off as ^ r'^\ causing £a to in- 
crease in magnitude, until the point at which the polarization 
effectively freezes. The distance from the star at which the 
adiabatic approximation breaks down is usually called the po- 
larization limiting radius (Heyl et al. 2003), and is given by 
the radius at which the relation B/\k- VB\ ~ Ia is satisfied 
(see Heyl et al. 2003; Lai & Ho 2003b for sUghtly different 
criteria). The coupling of the modes is gradual, however, and 
a characterization of the observed polarization signal requires 
solution of equations (21)-(22) over this coupUng region. 



It is instructive to quantify the size of the regions where 
resonant scattering and mode coupling are relevant, as this 
determines the simplifications that are possible in the full 
problem. The surface inside which resonant cyclotron scat- 
tering can occur is the escape radius, given by the condition 
(uic/uj)^ + = 1 (Fernandez & Thompson 2007). Outside of 
this surface, there is no particle velocity such that the res- 
onance condition (eq. [8]) can be satisfied. For the twisted 
dipole geometry, it is given by 



■ '■esc(^,'^)^)-Bpole, A^jv-S)^NS) 
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where Wc.poie = efipoie/ (W(,c), and 2^ = [F^(cos0)]'''^ is the an- 
gular dependence of the magnetic field strength (c.f. eq. [1]). 
The polarization limiting radius is set by the equality between 
£a and the length scale over which the field strength changes, 

eB=B/\k-VB\. Writing 4 = r/C(6', ^,A(/)n-s), where C is a 
dimensionless function of order unity, we find 



'■pi - '■pl(^)'^)^)5pole) A0N_s,/?Ns) 

2 e2 



= Rns 



146 



^pole N 2x^ f -^NS^ 

V307ry V^QED^ ^ ''^C V c 



Rns, 



-\ l/(i+2p) 

(27) 
(28) 



A<^=1 



where = /?ns / 1 0* cm. 

Figure 2 shows r^sc and r^i as a function of photon direc- 
tion angles, for 1 keV photons at the magnetic equator, with 
a twist A(/)NS = 1 and a standard choice for other parame- 
ters. Aside from narrow regions where the surfaces become 
singular (|/i| 1 for resc and k ■ VB for rpi), they are 
well separated by an order of magnitude in radius. Photons 
in the Rayleigh-Jeans tail of the thermal spectrum could in 
principle violate this separation, but they contribute little to 
the observed signal. One can thus assume that (1) scattering 
is largely decoupled from polarization evolution, taking place 
mostly within the adiabatic propagation region, and (2) the 
polarization signal is set at distances > 100/?ns from the cen- 
ter, where the effects of general relativity are weak. 



2.2. Seed Photon Distribution 

Two physical processes have been proposed to explain the 
thermal component of the quiescent soft X-ray spectrum. The 
first is heat generated in the core by ambipolar diffusion 
(Thompson & Duncan 1996). These photons are expected 
to emerge in E-mode polarization, as their free-free opacity is 
suppressed relative to the O-mode (e.g., Lai & Ho 2002). In- 
cluding resonant mode conversion at magnetar field strengths 
(> 10^"* G at the surface) does not alter this dominance of the 
E-mode, but instead modifies the spectral shape (e.g., Lai & 
Ho 2003a). For the purposes of exploring the effects of mag- 
netospheric scattering on the surface polarization, we adopt a 
simple blackbody spectrum with a temperature such that its 
value measured at infinity is = 0.4 keV. This case is as- 
signed a 100% initial polarization in the E-mode. 
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Fig. 2. — Escape radius rcsc (eq. [26], red surface) and polarization lim- 
iting radius r^i (eq. [27], blue surface) for 1 keV photons at the magnetic 

equator, as a function of photon direction 1(9^,4)]^), with 9]^ and (fyj^ the polar 
and azimuthal angle of k relative to the magnetic axis, respectively. Other 
parameters are iJpde = lO''* G, A(/>n-s = 1 rad, and i?Ns = 10* cm. The spikes 
in the surfaces arise at points where the expressions become singular, namely 
^1 — > 1 for /'esc and k ■ V6 for r^\. 



The second source of thermal photons is the deposition of 
energy at the stellar surface by returning currents in a layer 
that is optically thick to free-free absorption (Thompson et al. 
2002; Beloborodov & Thompson 2007). In this case, the 
emission is expected to emerge predominantly in O-mode po- 
larization. Since the details of this mechanism have not been 
thoroughly explored yet, we account for this possibility by 
simply changing the outgoing polarization distribution from 
E- to O-mode. Intermediate mixtures (such as a completely 
unpolarized case) can be obtained by linear superposition. 

Regardless of the mechanism responsible for the seed pho- 
tons, emission is expected to be anisotropic over the stellar 
surface. For the deep cooling case, even a constant temper- 
ature atmosphere has an emerging radiation intensity that is 
beamed relative to the magnetic field direction, with an an- 
gular distribution that depends on photon energy and mag- 
netic field strength (van Adelsberg & Lai 2006). Additional 
anisotropics are expected due to the suppression of thermal 
conduction across magnetic field lines. If a strong internal 
toroidal field is present, very narrow polar caps are expected 
(Perez-Azorm et al. 2006; Geppert et al. 2006). If the ther- 
mal emission is powered by returning currents, localized hot 
spots of the size of the current carrying region would be natu- 
ral. We explore the effects of anisotropics in the surface tem- 
perature by running additional models restricting emission to 
polar caps with a size of 5 degrees. We do not attempt here 
to include the effect of beaming of radiation around the local 
magnetic field direction. It is worth keeping in mind, though, 
that general relativistic deflection of photon trajectories sig- 
nificantly smears out emission from a polar cap, enlarging its 
apparent angular size (§2.3). 

2.3. Spacetime Properties 

Magnetars are slow rotators, with periods falling in the 
range 2- 10 s (e.g.. Woods & Thompson 2006). At the surface 
of the star, corrections to the metric due to rotation are of order 

-8 p2 



mass, and angular velocity, respectively, and Psec = P/l &■ 
Hence, using a Scharzschild spacetime is a very good approx- 
imation. 

In §2.1.3 we showed that the polarization signal is deter- 
mined at radii > IOOTJns, where the effects of general rela- 
tivity on photon propagation are weak. However, a significant 
fraction of the seed photons leaves the magnetosphere without 
scattering, thus inclusion of light-bending is needed to prop- 
erly account for the fraction of the stellar surface contributing 
to the unscattered polarization signal. This is particularly im- 
portant when including emission from hot spots, as bending 
angles are as large as 30 degrees for typical neutron star pa- 
rameters (e.g., Beloborodov 2002). 

We adopt the approach of Heyl et al. (2003) for inclusion 
of general relativity in our polarization calculation. That is, 
(1) modifications to the magnetic field are ignored, as close 
to the star this is of weak importance for both scattering and 
polarization^; (2) the polarization plane is parallel transported 
along photon geodesies; and (3) we ignore general relativistic 
corrections to the polarization transfer equations (21-22) other 
than the modification of the polarization plane. 

3. NUMERICAL IMPLEMENTATION 

Our Monte Carlo code has been extensively described in 
Fernandez & Thompson (2007). We describe here the modi- 
fications introduced in the present study. 

3.1. Propagation 

Integration of photon trajectories in a Scharzschild metric 
is a standard textbook calculation (e.g., Schutz 2009). How- 
ever, for practical purposes, we recast the equations in a form 
more amenable for the Monte Carlo algorithm. The calcula- 
tion proceeds as follows. For each spatial step of length A£, 
the coordinate system is rotated to a frame where x= r and 
z = r X k, e.g., the usual 6 = tt/2 plane (see Appendix A for 
the explicit transformation). We then integrate, as a function 
of azimuthal angle, the geodesies for the radial coordinate r, 
the components of ^ in a non-coordinate orthonormal polar 
basis, k = k''r + l 



(JNs/[M^sRNSc]f ^ (RNsi^/c) < 4 X 10-«/?^p-^„ where 
Jns, A/ns, and 57 = In/P are the stellar angular momentum. 



t>, where r = y/{l-rs/r)er and (f) = e^/r 
(e.g., Schutz 2009), and the line segment traversed Ai. Here 
and throughout the paper, = 2GMns/c^. The ordinary dif- 
ferential equation system is 



(29) 
(30) 
(31) 
(32) 



We have tested our implementation by verifying agreement 
with the asymptotic expansion of Beloborodov (2002), which 
is accurate to about 1 part in lO"-'. 

' It would affect scattering for either very high energy photons, very en- 
ergetic particle distributions, or ion scattering. In the case of polarization 
transfer, propagation close to the star is well inside the adiabatic limit; see 
§2.1.3. 
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Within the code, the problem is of boundary value type, as 
we require the line segment traversed As to equal the exter- 
nally imposed step size A£. This is needed so that a back- 
wards and forward step sequence (e.g., for spatial deriva- 
tives of the magnetic field) returns to the original position. 
In practice, we have found that a simple 3-step iterative pro- 
cess achieves this equality to double precision in most cases. 
The initial angular interval for the Runge-Kutta integrator is 
estimated as A0i = Alj ^ + (dr/ d^)^. Subsequent correc- 
tions are then A0i+i = (A£- Ai,) / ^r? + (dr/d0)?, where / is 
the order of the iteration. Despite these simplifications, com- 
bined solution of equations (21)-(22) and (29)-(32) becomes 
computationally expensive when > 10^ photons are followed. 
We thus cut off light bending when r^/r < 2 x 10"^, i.e., when 
general relativity modifies the trajectory to less than 1 %, and 
impose planar propagation outside that radius. 

3.2. Scattering 

Aside from taking /3 out of the square brackets in equa- 
tion (12), the scattering algorithm remains the same as in Fer- 
nandez & Thompson (2007), including the adaptive step size 
algorithm in momentum space. 

The overlap of the photon polarization with a circularly po- 
larized wave |e±p, which enters the optical depth (eq. [12]), 
needs to be modified for the possibility of non-adiabaticity. 



1 



A - ( flrX- \/ 1- firZ±iy 

2 



= -[M,Re(A,)Tlm(A,)] 
+i[/i,Im(A,)±Re(Av)]', 



(33) 



where k = z, and the magnetic field lies in the x-z plane. Thus, 
when modes are coupled, the sign of the charge needs to be 
taken into account. 



3.3. Polarization Evolution 

Since £a oc B~^, this length scale can become very small 
at radii r ^ rpi, hence integration of equations (21)-(22) be- 
comes impractical in regions where adiabatic propagation is a 
good approximation. An arbitrary choice has to be made for 
the transition surface where the ordinary differential equation 
solver is turned on. Once this choice is made, we set the initial 
conditions for the integration such that the absolute value of 
the components of the electric field correspond to those of the 
normal mode the photon was in before reaching the transition 
surface. A single random complex phase is then assigned to 
all components. 

Figure 3 shows the real part of the electric field vector rel- 
ative to the magnetic axis for a test E-mode photon moving 
outwards along the x-a\is, in a magnetosphere with a twist 
A(/)Ns = 1- Integration of equations (21)-(22) begins at 



^couple 



^ = 10-1 



(34) 



corresponding to a radius r ~ 407?ns on the equatorial plane. 
As the photon travels out, the characteristic spatial oscillation 
frequency of A becomes longer, and the polarization reaches 
its asymptotic value at a distance ~ 200/?ns from the star In 




1000 



Fig. 3. — Real components of the electric field (solid lines) and normal 
mode eigenvectors (dashed-lines) for a test E-mode photon moving radially 
outward along the x-axis. Integration of equations (21)-(22) begins when 
'?coiiplc = 10 (sq- [34]). The photon stays in the same normal mode because 
the component of B transverse to the propagation does not change direction. 



this particular case, because the trajectory is radial, the com- 
ponent of the magnetic field perpendicular to the trajectory 
does not change direction, thus normal modes should not cou- 
ple (c.f. eq. [24]). Figure 3 shows that the amplitude of the 
oscillations in A are bracketed by the absolute value of the 
corresponding polarization eigenvector, hence our integration 
correctly captures adiabatic propagation. This purity of nor- 
mal mode state is achieved to double precision for radial prop- 
agation, as measured by the intensity in the orthogonal mode. 

Figure 4 shows the effect of changing the value of ?7coupie for 
a test E-mode photon moving along r = IST^ns-^ + A(i +}'+£), 
where A > is a parameter labeling the trajectory^. Because 
adiabatic propagation is captured with good accuracy, there is 
no significant difference in starting the integration as deep as 
?7coupie = 10"^, corresponding roughly to a factor two in dis- 
tance relative to ?7coupie = lO""*, showing that the final photon 
state is converged. 

As the photon moves out, the polarization tensors 
(eqns. [18]-[19]) smoothly approach unity, with a consequent 
freezing of the electric field vector Figure 4a shows that the 
absolute value of the spatial components of A relative to the 
magnetic axis cease to vary significantly at radii ^ 200/?ns- 
Nonetheless, a quantitative criterion is needed to establish 
freezing in practice. As a figure of merit, we measure the 
norm of the difference between the polarization vector before 
and after each spatial step 



|AA| = |A(r + Aa)-A(r) 



(35) 



and then proceed to multiply this quantity by the factor 
{r/ At), obtaining an estimate of the change in A over a scale 
height in magnetic field strength. The result is shown in Fig- 
ure 4b. The dashed line shows the point where \AA.\(r / At) = 
10"^ which we establish as our fiducial criterion for polariza- 
tion freezing. As shown by Figure 5, this rather conservative 
limit is located at distances < 1000/?ns from the star for most 
trajectories, safely below the light cylinder by at least an order 
of magnitude. It also lies below the point where the plasma 

We ignore Hght-bending in Figures 4 and 5. 
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Fig. 4. — Various quantities for a test E-mode photon moving along the line 
18RnsA'+ ACi+y+l) relative to the magnetic axis, where A > is a parameter 
labeling the trajectory. Different curves correspond to different values r^^oupic 
above which the polarization evolution equations (21)-(22) are solved. Panels 
show, as a function of radial distance from the star, (a) the absolute value 
of the electric field vector components relative to the magnetic axis, (b) the 
approximate change of the polarization vector over one magnetic field scale- 
height, (c) the projection of the electric field onto the 0-mode eigenvector, 
and (d) the normalized Stokes V parameter (eq. [39]). The dashed line shows 
the value of the electric field change at which we consider the polarization 
frozen (see text for details). 




10000 
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Fig. 5. — Norm of the change in the polarization vector over a scale height 
in magnetic field (eq. [35]) for vaiious trajectories, labeled by the direction 
of the wave vector k relative to the magnetic axis. The latter is parameterized 
in terms of the polar and azimuthal angles 9j- and respectively. Photons 
start from r = 40Sns^> and their wave vector points towards x (black), [x + 
y + ?] (red), z (green), and [-.v + |] (blue). The latter trajectory crosses the 
z-axis. The horizontal dashed line marks the point where we consider the 
polarization to be effectively frozen. 




J 2% 



contribution would become comparable to that of the magne- 
tized vacuum (eq. [17]). 

Figure 4(c)-(d) also illustrates how normal modes become 
coupled in a non-radial trajectory. Panel (c) shows the nor- 
malized intensity in the O-mode, \Ao\^ = |A • eoP with eo the 
O-mode eigenvector, while panel (d) shows the normalized 
Stokes V component (eq. [39]). IVIodes begin to couple within 
a factor 2 of the polarization limiting radius (eq. [27]), before 
the electric field components freeze. This particular trajectory 
asymptotes to a ^ 10% amplitude in the orthogonal mode, 
and a circular polarization component of the same magnitude, 
showing the importance of solving equations (21)-(22). 

As a further test of our implementation, we have calculated 
the polarization signal as a function of phase for two configu- 
rations with = (no scattering), light bending, and other 
parameters mirroring Figures 7 and 8 of Heyl et al. (2003)^. 
The latter are equivalent to a configuration with stellar ra- 
dius 7?NS = 10 km, a polar magnetic field Bpoie = 2 x 10'^ G, 

' The panel in the lower left could only be matched to the corresponding 
Figure in Heyl et al. (2003) by changing the orientation angles. 



Fig. 6. — Polarization fraction for monoenergetic photons and dipolar field 
(no scattering) as a function of phase for different relative orientations of 
magnetic axis, rotation axis, and line-of-sight (eqns. [42] and [43]). Param- 
eters are Rns = 10 km, Bpoie = 2 X lO'- G, and tiuj = 6.626 eV (black) and 
0.6626 keV (red). Results agree with Figures 7 and 8 of Heyl et al. (2003). 



and monoenergetic photons with energy huj = 6.626 eV and 
0.6626 keV, respectively. 

3.4. Postprocessing 

Once the photon polarization is frozen, we store its fre- 
quency, position, direction, and polarization vectors. In con- 
trast to Fernandez & Thompson (2007), we do not convolve 
a frequency shift distribution with a blackbody function, but 
instead directly draw photons from a Planck distribution. The 
convolution of a frequency independent shift distribution with 
a source function is incompatible with the fact that the polar- 
ization evolution equations depend explicitly on the photon 
frequency through A;o (eqns. [21]-[22]). 

In addition to the aforementioned variables, we keep track 
of the number of scatterings a photon has had before escaping. 
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TABLE 1 

Model Parameters 



Name" 


Dist.'' 


A(^N(; (rad) 


7max 
max 


Seed Pol. 




Fis. 


tl0g20c90Eu 


uni 


1 


2 


E-mode 


90 


7-20 


tl0g20c90Eb 


bid 


1 


2 


E-mode 


90 


12 


t08 g24c9()Eu 


uni 


0.8 


2.4 


E-mode 


90 


13,18 


tl2gl8c9()Eu 


uni 


1.2 


1.8 


E-mode 


90 


13,18 


tl4gl7c90Eu 


uni 


1.4 


1.65 


E-mode 


90 


13,18 


tl6gl6c90Eu 


uni 


1.6 


1.6 


E-mode 


90 


13,18 


tl0g20c90Ou 


uni 


1 


2 


0-mode 


90 


15,19 


tl0g20c5bEu 


uni 


1 


2 


E-mode 


5 (both) 


16,20 


tl0g20c5sEu 


uni 


1 


2 


E-mode 


5 (south) 


16,20 


tl0g20c5pEu 


uni 


1 


2 


E-mode 


5 (planar) 


16,20 



"We do not list here models with Ac^ = 0, or test calculations. 
''Directionality of the particle energy distribution (eq. [7]), with uni denoting uni- 
directional, and bid bidirectional, respectively. 



and construct separate histograms labeled by this parameter. 

The standard parameterization of polarization degrees of 
freedom is a set of four Stokes parameters (e.g., Rybicki & 
Lightman 2004). Choosing the photon direction as the z axis, 
there is an additional degree of freedom in their definition: the 
rotation angle around this axis. Choosing a coordinate system 
where the magnetic field is initially in the x-z plane, we can 
write 



I=\A^\^+\Ay\^ 

e=|A,|2-|A,|2 
f/ = 2Re(A,.A*) 
y = 2Im(A,A;), 



(36) 
(37) 
(38) 
(39) 



where the star stands for complex conjugation, and / measures 
the intensity, Q and U linear polarization, and V circular po- 
larization. Standard observables are the linear polarization 
fraction 



and the polarization angle 



1 fu 
Y = - arctan — 

2 \Q 



(40) 



(41) 



Rotation of the Stokes parameters around the z-axis is 
straightforward (e.g., Chandrasekhar 1960), and involves only 
Q and U. For phase dependent quantities, we perform a rota- 
tion such that the x-axis is in the plane of the line of sight and 
the rotation axis. 

Quantities as a function of phase are obtained by perform- 
ing suitable geometric transformations (Appendix A), which 
have as input parameters the two orientation angles 



cos9mt = M -fl 

COS^los = 0-0, 



(42) 
(43) 



where M, Cl, and O are the magnetic axis, rotation axis, and 
line of sight direction, respectively^'. 

We focus on linear over full polarization fraction because photoelectric 
polarimeters will not be sensitive to circular polarization (e.g., Swank et al. 
2010). 

" These two angles are denoted respectively by ^ and x in Nobili et al. 
(2008). 



3.5. Model Parameters 

Table 1 shows the sequence of models investigated in this 
paper. All of them cover the energy range 0.04-40 keV. We 
run two versions of each. One with 40 bins per decade in 
energy and 64 bins in cos 0k, optimized for angle-averaged 
studies as a function of energy. We also compute models 
with 10 bins per decade in energy and 256 bins in angle, to 
compute frequency integrated lightcurves. Because the mag- 
netic field is axisymmetric, we average our results in cf)/^. The 
number of seed photons is such that there are 10^ per energy- 
angle bin. The stellar radius is assumed to be /?ns = 10 km, 
the Schwarzschild radius /?ns /''s = 3, the magnetic field at the 
poles Bpoie = lO'"* G, and the temperature of the seed black- 
body measured at infinity kT^o = (1 - '"s/^ns)'''^ kT = 0.4 keV. 
Other parameter combinations for each model are Usted in 
Table 1. Model names refer to their twist angle (t), maxi- 
mum Lorentz factor (g), polar cap radius (c), seed polarization 
(E/O), and directionality (u/b). 

4. RESULTS 

Since one of our primary goals is to provide predictions for 
future X-ray polarimeters like GEMS, we focus on the behav- 
ior of the linear polarization fraction (eq. [40]) and polariza- 
tion angle (eq. |41|) in relation to the total intensity, as dif- 
ferent system parameters are changed. Phase-integrated and 
phase-dependent observables provide complementary infor- 
mation about the system, and we consider each class of mea- 
surement in turn. Due to rotation, both sets of observables 
depend on the relative orientation of rotation axis, magnetic 
axis, and line-of-sight, as parameterized by the angles O^ot and 
^^los (eqns. [42] and [43], respectively). Therefore, we begin 
by discussing the emission in the corotating frame, where the 
physics giving rise to the polarization signal is most easily 
understood. 

4.1. Polarized Emission in the Corotating Frame 

We first focus on the basics of the angular and frequency 
distribution of the polarized emission in the corotating frame. 
Figure 7 shows two-dimensional histograms of intensity, po- 
larization fraction, and polarization angle measured at infin- 
ity, as a function of photon energy and magnetic colatitude, 
for model tl0g20c90Eu. For illustration, we have separated 
out in different columns the contributions from scattered and 
unscattered photons, plus their sum. 
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Fig. 7. — Observables measured at infinity for model tl0g20c90Eu, as a function of photon energy and magnetic colatitude. Each panel shows, from top 
to bottom, number of photons collected, linear polarization fraction (eq. [40]), and polarization angle (eq. [41]), respectively. Left, middle, and right panels 
correspond to the total, unscattered, and scattered contributions, respectively. Regions of parameter space with no counts are shown in black. 



Panels (a)-(c) of Figure 7 show the number of photons col- 
lected at infinity. The unscattered distribution has an energy 
dependence close to that of a blackbody'^ and an angular de- 
pendence given by /unscatt = Ioe~'^^^\ where Iq is the surface 
intensity and r ~ sin^ 9A(j)//3 the total optical depth (the inte- 
gral along the line-of-sight of eq. [12]). The exact dependence 
of the latter on magnetic colatitude and the unidirectional par- 
ticle velocity distribution introduce a slight north-south asym- 
metry. The scattered component, on the other hand, peaks at 
the southern hemisphere. This results from the fact that pho- 
tons with a large positive energy shift move out of the magne- 
tosphere when they are generated close to the magnetic equa- 
tor, where the magnetic field points downwards. Upscattered 
photons moving northward are re-entering the magnetosphere 
(see below). At fixed magnetic colatitude, the frequency de- 
pendence is smooth curve and extends to both sides of the 
original blackbody peak. 

Panels (d)-(f) of Figure 7 show the corresponding polar- 
ization fractions. The unscattered component reflects the po- 

The optical depth (eq. [12]) depends implicitly on frequency through 
the resonant momenta (7/?)^. 



larization distribution of seed photons, in this case 100% E- 
mode. The scattered component has a characteristic polar- 
ization fraction of ^ 30%, which arises mostly from mode 
exchange (see below). The total polarization fraction peaks 
at the north pole, because more unscattered photons leave in 
that direction, and they are more polarized. This is the con- 
sequence of having a uni-directional particle energy distribu- 
tion, which causes a break of the north-south symmetry. Thus, 
the polarization fraction carries a signature in both energy and 
angle. 

From equation (16), one can see that photons are more 
likely to be in the E-mode state immediately after scattering, 
for all < 1. This, coupled with nearly adiabatic propaga- 
tion, means that E-mode polarization dominates the scattered 
emission for all values of the energy and magnetic colatitude. 
If photons were to scatter only once, the net polarization frac- 
tion would be 



1-(K) 



l + (/i;)2 



(44) 



where the weight [l+(/ir)^] comes from the total differen- 
tial cross-section (eq. [13]). The characteristic polarization 
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Fig. 8. — Top: Polai'ization fraction in tlie first scattering order for a run 
witli tile same parameters as modei tlO g20 c90 Eu, but forcing pliotons to es- 
cape after one scattering. Middle and Bottom: Fraction of re-scattered E- and 
0-mode pliotons, respectively, obtained by subtracting the normal mode in- 
tensities of model tlO g20 c90 Eu (which undergoes multiple scattering) from 
the corresponding values from the ran shown in the top panel (with a single 
scattering). As E-mode photons have a larger cross-section, they undergo rel- 
atively more re-scatterings than O-mode photons and hence the polarization 
fraction decreases from the upper limit given by equation (44). 



fraction of ^ 30% observed at infinity arises from multiple 
scattering. In particular, E-mode photons have a larger cross- 
section than O-mode photons at all angles because of the over- 
lap \e±\- (eqns. [14] and [33]), and hence they are more likely 
to undergo further scatterings, decreasing the overall polariza- 
tion fraction. 

To show this effect. Figure 8(a) displays the polarization 
fraction of the scattered component in a run with the same 
parameters as model tl0g20c90Eu (Figure 7) but this time 
forcing photons to escape after their first scattering. The aver- 
age polarization fraction is consistent with what equation (44) 
predicts. One can take the difference between the normal 
mode intensities in this run and those from the first scatter- 
ing order in model tl0g20c90Eu to find the fraction of re- 
scattered photons. The resulting numbers for E- and O-mode 
are shown in panels (b) and (c) of Figure 8. As expected, 
E-mode photons undergo relatively more re-scatterings than 
O-mode photons, decreasing the overall polarization frac- 
tion of the first scattering order while still dominating the 
output. The bulk of re-scattering affects photons that move 
towards the north magnetic pole and have a high energy. 
The maximal frequency boost imparted by equation (8) is 
(1 +/3/iin)/(l -/3/iout), with /^in > and ^i^m > the ingoing 
and outgoing cosine between photon direction and magnetic 
field (/3 > in the uni-directional case). Hence these pho- 
tons correspond to those who first scattered closer to the south 
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Fig. 9. — Circular polarization fraction measured at infinity (normalized 
Stokes V parameter, eq. [39]) for model tlO g20 c9() Eu, as a function of pho- 
ton energy and magnetic colatitude. 



magnetic axis, where the magnetic field points mostly north- 
ward, and whose value of ^ changed from negative to positive. 
After their first scattering, these photons move inwards, being 
most likely to scatter again. An analogous argument shows 
why high-energy photons going southward after scattering are 
the most likely to escape. 

Panels (g)-(i) of Figure 7 show the polarization angle. It is 
immediately evident that this quantity is mostly independent 
of energy for this particular case, varying only with magnetic 
colatitude, aside from a small region of scattered emission 
near the south pole. Absent twist, the magnetic field is dipolar 
and parallel to the magnetic axis along the line-of-sight This 
would yield a polarization angle of 90° for the E-mode and 
0° for O-mode. A positive A(/)ns introduces a dependence on 
magnetic colatitude, increasing these values away from the 
poles. The functional form is equal, within statistical errors, 
to the angle that the magnetic field makes with the magnetic 
axis. For model tl0g20c90Eu and measured relative to the 
magnetic axis, the polarization angle is x = arctan(B0/B0) + 
tt/2. 

The weak dependence on energy results from E-mode po- 
larization dominating the scattered component at most angles 
and energies, combined with the assumption of E-mode seed 
photons in the tl0g20c90Eu model. This dominance of the 
E-mode in the scattered emission persists even if 100% O- 
mode seed photons are used. The only difference is a 90° 
jump as a function of energy when transitioning to the unscat- 
tered component. The polarization vector follows the local 
magnetic field direction adiabatically until the coupling re- 
gion at r ^ Tpi (§2.1.3). Hence, the polarization angle of the 
scattered component at most latitudes will trace the magnetic 
field direction at r ^ rpi ^ IOOT^ns (eq- [27]). 

Figure 9 shows the degree of circular polarization for model 
tl0g20c90Eu as a function of energy and magnetic colat- 
itude. It is produced almost exclusively by the scattered 
component, and its value deviates from zero by a few per- 
cent throughout the energy-colatitude plane, except in regions 
close to the magnetic poles at high energies. The maximum 
absolute value of about 8% is reached near the south pole, 
towards which most scattered photons with a positive energy 
shift escape. 

In principle, any photon moving non-radially will develop 
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Fig. 10. — Intensity (top), polaiization fraction (middle), and polarization 
angle (bottom) as a function of rotational phase <p^„t in two energy bands 
(2-4 keV and 4-12 keV) for model tl0g20c90Eu and various viewing ge- 
ometries. The curves coiTespond to pairs of angles (9rot, ^*ios) of (45°, 70°: 
solid/black), (70°, 45°: dotted/red), (60°, 70°: dashed/blue),' and (90°, 90°: 
dot-dashed/green). For illustration, two full periods are shown. 



a circular polarization component as it traverses the magneto- 
sphere (e.g., Figure 4[c]). To see this, one can write an equa- 
tion for the rate of change of V along the photon trajectory, by 
combining equations (21), (22), and (39). After some algebra, 
one finds that to leading order in 6 = Q;em(fi/fiQED)^/(457r) one 
has 

dV T n 

— = {q + m)sm^ 9kBCOs2(j)kBU + 0(S^), (45) 

oz 

where 6*^^ and are respectively the polar and azimuthal 
angle of the magnetic field in a coordinate system where k = z 
and B is in the x-z plane. If the polarization vector is in a 
normal mode state, then U = and no circular component is 
generated. Similarly, dV/dz = if photons propagate parallel 
to the magnetic field. 

The fact that the unscattered component has V almost iden- 
tically zero everywhere is due to the fact that a given non- 
radial trajectory originating from the star always has a coun- 
terpart for which the magnetic field rotates in the other di- 
rection along the trajectory, cancelling out the net contribu- 
tion to V if they have the same intensity. The scattered com- 
ponent, however, is not isotropic and thus such cancellation 
does not occur, leaving a residual circular polarization sig- 
nal. If plasma polarization were to become dominant at some 
radius comparable to rpi (eq. [27]), then the circular polariza- 
tion fraction could be much higher and substantially modify 
the outgoing linear polarization signal. 
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Fig. 1 1 . — Magnetic colatitude swept by the line-of-sight as a function of 
rotational phase for vaiious viewing geometries. The color and line styles 
have the same meanings as in Figure 10. 



4.2. Phase-resolved Emission 

Ideally, one would like to observe the system from all an- 
gles in the corotating frame and generate observed maps as 
a function of magnetic colatitude for comparison with plots 
like Figure 7. In practice, we have a fixed line of sight for 
any single source, so that the viewing angle corresponds to a 
time varying magnetic colatitude as a function of rotational 
phase (proi- Due to this relatively direct relationship with the 
corotating frame observables, we find it pedagogically advan- 
tageous to discuss phase-resolved observables first. This rela- 
tionship also means that polarization measurements as a func- 
tion of rotational phase yield the most detailed information 
about the system. Hence, they are the most useful observa- 
tional constraints for disentangling geometrical effects from 
the physics of the surface and magnetosphere. However, they 
may necessitate much longer integration times to achieve use- 
ful signal-to-noise. We discuss prospects for obtaining such 
phase-resolved measurements with GEMS in §5. 

Figure 10 shows several light curves corresponding the to 
the tlO g20 c90 Eu model for different lines-of-sight and mag- 
netic axis orientations. The phase is defined so that the north 
and south magnetic poles cross the line-of-sight at 0rot = 
and TT, respectively. Two full periods are shown for illustra- 
tion. The intensity, polarization fraction, and polarization an- 
gle are computed after integrating the emission over a soft 
energy band (2-4 keV) that is dominated by direct emission, 
and a hard band (4-12 keV) which is mostly scattered flux. 
The normalization is chosen so that the average over phase is 
unity. 

The variation of the total intensity and polarization fraction 
follows directly from the two dimensional map in Figure 7, 
as there is a simple mapping between rotational phase and 
magnetic colatitude (eq. [A7]). Figure 1 1 shows this mapping 
for the same set of orientations (6'rot, 6*108) that are plotted in 
Figure 10. Orientation effects depend mostly on the maxi- 
mum (most negative cos 9) observable colatitude, because the 
polarization fraction tends to be a minimum and the total in- 
tensity a maximum near the south magnetic pole. The map- 
ping is invariant under the exchange 6'rot ^ 6*105, so the solid 
black (6l„t = 45°, 6iio,s = 70°) and dotted red curves (e.^t = 70°, 
6*105 =45°) are identical. 

The orthogonal rotator viewed at 6*105 = 90° gives the largest 
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Fig. 12. — Intensity (top) and polarization fraction (bottom) as a function 
of phase in two energy bands (2-4 keV and 4-12 keV). The curves correspond 
to models tl0g20c90Eu (black solid) and tl0g20c90Eb (red dotted). The 
former assumes that the particle distribution is unidirectional, while that of 
the latter is bidirectional. Magnetic and viewing orientations corresponding 
to 0rot = 45° and 9uk = 70° are assumed. 



modulation. In this case the observer sees one of the mag- 
netic poles directly at (pj-ot equal to integer multiples of tt. The 
optical depth to scattering is lowest along the poles, leading 
to dips in the total intensity at these phases. The polarization 
fraction is also low since the photon wave vector lies parallel 
to the magnetic field and there is no preferred direction for the 
polarization. 

The gross, qualitative features of the total intensity and po- 
larization fraction phase curves are largely insensitive to the 
orientation. There is a strong peak near 0iot = tt and generally 
a weaker one near 0iot = 0, when the south and north mag- 
netic poles (respectively) cross the line-of-sight. The contrast 
is particularly clear in the hard energy band which is domi- 
nated by the scattered emission. However, in the soft band, 
which is dominated by the unscattered component the polar- 
ization fraction peaks near ^i-d = and has a minimum near 

01-ot = TT. 

The half period phase shift between the polarization frac- 
tion and the total intensity results from the assumption of a 
unidirectional photon distribution. Figure 12 shows a compar- 
ison of the tl0g20c90Eu and tl0g20c90Eb models. These 
differ only in their treatment of the particle distribution, with 
the latter having a bidirectional charge distribution. The 
largest contrast arises when comparing the intensity in the 
hard band with the polarization fraction in the soft band. The 
intensity and polarization fraction of the bidirectional model 
is completely symmetric in magnetic colatitude, so the differ- 
ences between (p^oi = and vr are purely orientation effects in 
this case'^. In contrast to the unidirectional model, the peaks 
in intensity and polarization fraction coincide. Thus, one can 
distinguish whether the magnetosphere is filled with electron- 
positron or electron-ion pairs because this difference in parti- 
cle distribution manifests itself as a break in the north-south 
symmetry for the unidirectional model. 

As mentioned in §4.1, the polarization angle is determined 
primarily by the projection of the magnetic field at r ^ r^i 
onto the line-of-sight. In phase resolved plots, we have cho- 
sen the zero point for this angle to be the rotation axis, so 

The observable colatitude range is [25°, 115°] 
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Fig. 13. — Intensity (top), polarization fraction (middle), and polariza- 
tion angle (bottom) as a function of phase in two energy bands (2-4 keV 
and 4-12 keV) for different twist angles. The curves correspond to mod- 
els tl0g20c90Eu (solid, black), t08g20c90Eu (dotted, red), tl2g20c90Eu 
(blue, dashed) tl4g20c90Eu (green, dot-dashed), and tl6g20c90Eu (vio- 
let, triple-dot-dashed). Magnetic and viewing orientations corresponding to 
^'roi =45° and Oi^s = 70° are assumed. For comparison, the long-dashed black 
line shows the polarization angle as a function of phase for a run with no twist 
and thus a purely dipolar field. 



the absolute value will differ from that shown in Figure 7 
by a geometric transformation involving 0,ot, ^los, and the ro- 
tational phase 4>rot (eqns. [A7]-[A10]). This means that the 
phase curve of polarization angle is not affected by the in- 
variance to O^oi ^ dios exchange of the phase-colatitude map- 
ping, as there is an additional rotation of the polarization vec- 
tor as the magnetic axis rotates around the line-of-sight. This 
is seen most clearly by comparing the (6',ot = 45°,(?io,s = 70°) 
and (6'i-ot = 70°, 6*105 = 45°) curves in Figure 10. Although the 
intensity and polarization fraction phase curves are identical, 
the variation with polarization angle differs greatly between 
the two orientations. In general, for 6'rot < 6*105, the polarization 
oscillates about a mean value, while for 6'rot > 6*10,5 it sweeps 
through the full allowed range of 180°. This has important 
implications for the phase-averaged polarization, which we 
discuss below. 

Figure 13 shows the variation of the phase-resolved observ- 
ables as twist angle A0ns and maximum Lorentz factor 7inax 
of the momentum distribution (eq. [7]) are varied. We are 
primarily interested in the variation with A(/)ns, but we vary 
7max simultaneously (see Table 1) to keep the power law slope 
of the spectrum as uniform as possible. As twist angle in- 
creases, the intensity increases for all energies above ^ 2 keV 
(c.f. §4.3) as the fraction of scattered photons increases. Thus 
the peak in the intensity near (j>rot = ti" becomes stronger in both 
bands. 
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Fig. 14. — Polarization angle as a function of phase for a model with no 
twist (no scattering) and 100% E-mode seed photons for the same orientations 
as in Figure 10. Results for the soft and hard bands are identical due to lack of 
scattering, hence only one panel is shown. Note that the curves are symmetric 
every 1/2 phase due to the dipolar geometry of the zero-twist field. A net 
twist angle of the field lines introduces an asymmetry for most orientations. 



The polarization fraction shows the opposite trend with 
twist angle, however As A(/)ns increases, the polarization 
fraction gets smaller at every phase in both the soft and hard 
bands. For the soft band, the primary explanation is that scat- 
tered photons have lower polarization than unscattered, so that 
a larger fraction of scattered photons yields a lower overall po- 
larization fraction. As with the intensity variation, the effect 
is again largest at (j)rot = tt, when the south pole crosses the 
line-of-sight. This result is qualitatively similar to that found 
by Nobili et al. (2008) using adiabatic propagation and taking 
the relative difference between the normal mode intensities 
to obtain polarization fractions. Nobili et al. (2008) averaged 
over magnetic colatitude and energy to obtain their result, and 
offered a similar interpretation. 

However, it is notable that the drop in polarization fraction 
occurs even in the hard band. Unlike the soft band, it con- 
tains very little highly-polarized, unscattered emission, so the 
above explanation does not apply. Instead, the drop in po- 
larization with increasing twist angle is due to the increased 
number of scatterings per photon, which remove relatively 
more E-mode than O-mode photons from the escaping beam 
along a given direction (see also §4.1). 

For the orientation of Figure 13, the polarization angle os- 
cillates about a mean value near 130°, with an overall shape 
that depends weakly on twist angle but which becomes more 
asymmetric as this parameter is increased. For reference, the 
black dashed curves in the panel corresponding to the soft 
band shows the variation in the polarization angle with phase 
for a model with Ai^ns = 0. In this case the oscillation is sym- 
metric about 90° because the field is dipolar Hence, asym- 
metry in the polarization angle is a direct indication of the 
presence of a twist angle at r ^ rpi, and the magnitude of the 
asymmetry constrains the value of A^ns- 

To illustrate this effect in more detail. Figure 14 shows the 
polarization angle that is obtained when the twist angle is set 
to zero, for the same set of orientations shown in Figure 10 
(only one panel is shown, as the lack of scattering makes the 
results identical in both bands). Except for the orthogonal 
rotator, all curves in Figure 10 break the symmetry at 1/2 
phase relative to Figure 14, showing how the presence of a 
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Fig. 15. — Intensity (top), polarization fraction (middle), and polarization 
angle (bottom) as a function of phase in two energy bands (2-4 keV and 4- 
12 keV) as the seed photon polarization is varied. The curves correspond to 
models with seed distributions that are 100% E-mode (solid, black); 100% 
O-mode (dashed, red); 50% E-mode + 50% O-mode (blue, dashed); 25% E- 
mode + 75% O-mode (green, dot-dashed); and 75% E-mode + 25% O-mode 
(purple, triple-dot-dashed). Magnetic and viewing orientations correspond- 
ing to 9rot = 45° and 6ios = 70° are assumed. 



net twist angle in the magnetosphere becomes directly mea- 
surable through the polarization angle. We emphasize that this 
effect is a more general consequence of the properties of QED 
eigenmodes (Heyl & Shaviv 2000), and only indirectly due to 
magnetospheric scattering in that a definite polarization mode 
becomes dominant. It would be worth exploring under what 
circumstances it is possible to uniquely reconstruct the mag- 
netic field geometry from a given polarization angle curve; 
this work however lies beyond the scope of this paper. 

We now consider the dependence of observables on the seed 
photon polarization, which so far has been 100% E-mode. 
Figure 15 shows the variation of the phase-resolved quantities 
as the polarization of the seed photon distribution varies. Only 
two models with 100% E-mode and O-mode seed polarization 
are directly computed (tlO g20 c90 Eu and tlO g20 c90 Ou, re- 
spectively). The linearity of the Stokes parameters allows us 
to combine these runs to get any partial polarization. Here we 
consider an unpolarized (50% E-mode, 50% O-mode) seed 
distribution and a two partially polarized (75% E-mode and 
25% O-mode, 25% E-mode and 75% O-mode ) to illusti-ate 
the trends. The soft band intensity changes because the op- 
tical depth has a different colatitude variation in the pure O- 
mode model (tl0g20c90Ou), due to the angular dependence 
in the overlap function |e±p (eqns. [14] and [33]). The hard 
band has a slightly weaker peak for the pure O-mode model 
due to the reduced scattering fraction. The unpolarized and 
partially polarized models lie between the pure polarization 
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Fig. 16. — Intensity (top) and polarization fraction (bottom) as a function 
of phase in two energy bands (2-4 keV and 4-12 keV) for various assump- 
tions about NS surface emission and photon propagation. The curves corre- 
spond to full surface emission with light bending included (tl0g20c90Eu, 
solid black), 5° polar cap on both magnetic poles with light bending 
(tlOg20c5bEu, dotted red), 5° polai' cap with light bending but from the 
south pole only (tlO g20 c5s Eu, dashed blue), and 5° polar caps from both 
poles and planar propagation (tlO g20c5p Eu, dot-dashed green). Magnetic 
and viewing orientations con'esponding to Smt = 45° and 6\a^ = 70° are as- 
sumed. 
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Fig. 17. — Phase-averaged intensity (top), polarization fraction (middle), 
and polarization angle (bottom) as a function of photon energy for model 
tlOg20c90Eu and various viewing geometries. The curves correspond to 
the same orientation angles as in Figure 10. 



models by construction. 

The polarization fraction shows a somewhat stronger de- 
pendence on the seed photon distribution. For the pure O- 
mode model, the soft band polarization fraction is much 
weaker that the pure E-mode case at ^rot = tt. This is be- 
cause the transition from dominance of unscattered O-mode 
seeds to E-mode dominated scattered emission occurs within 
this band. As expected, the unpolarized and partially polar- 
ized models yield lower polarization fractions in the soft band, 
where unscattered photons dominate. The unpolarized model 
has a non-zero polarization fraction solely because of the non- 
negligible scattered emission towards lower energies. In the 
hard band, the scattered emission is E-mode dominated for all 
models, and the polarization fraction generally increases as 
the E-mode fraction of the seed photon increases. 

In the soft band, which is dominated by unscattered pho- 
tons, the phase curves of polarization angle for the models 
with pure E-mode and pure O-mode seeds differ by exactly 
90°. The partially polarized models follow the curve of the 
dominant polarization, while the unpolarized model switches 
between the two. In the hard band, the pure E-mode and pure 
O-mode models yield identical polarization angles at some 
phases and differ by 90° at others. These discontinuous jumps 
occur in the 100% O-mode seeds case due to a colatitude- 
dependent transition from O-mode seeds to E-mode domi- 
nated scattered photons. The partially polarized and unpo- 
larized models follow the E-mode curve. 

Thus far all the models have assumed uniform emission 
of seed photons from the NS surface. Figure 16 illustrates 
the effects of confining emission to polar caps of 5° radius 
around the magnetic poles. In this case we only show the 
intensity and polarization fraction, as the polarization angle 
curves are identical. The solid (black) curves correspond to 
emission from the full surface (tl0g20c90Eu) while the red- 
dotted curves refer to a model with polar caps around both 
magnetic poles (tl0g20c5bu). For the latter, there are mod- 
est enhancements and reductions of the intensity at (f)i-ot = 



and vr, respectively. This results primarily from the magnetic 
colatitude dependence of the scattering cross section, which 
is lowest near the poles. 

Although there are these modest differences, the phase 
curves are remarkably similar for the two models, given the 
large differences in the fraction of the star that emits. Even 
though the polar cap emission is confined to a very small 
solid angle, light bending causes the NS to appear as a more 
isotropic emitter for observers several stellar radii from the 
star where much of the scattering occurs. For comparison, the 
green (dot-dashed) curves show a model with planar propa- 
gation and no light bending (tl0g20c5pu). The differences 
with the isotropic emitter are much larger and there is much 
clearer signature in the polarization fraction. 

Finally, we consider a case in which only the south mag- 
netic pole radiates, including light-bending (tl0g20c5su). 
This corresponds to the blue (dashed) curve in Figure 16. 
Since we have chosen an orientation in which the observer pri- 
marily sees the northern hemisphere, the lightcurve peaks at 
0rot = TT, when the south pole crosses the line-of-sight. When 
the north pole crosses our line-of-sight at (p^oi = 0, there are 
very few unscattered seed photons, resulting in a low polar- 
ization fraction. 

4.3. Phase-averaged Emission 

Averaging quantities the over stellar rotation yields higher 
signal-to-noise. Hence, phase-averaged measurements are 
likely to be the first (or perhaps only) constraints obtained 
when polarimeters become operational. 

Figure 17 shows the tl0g20c90Eu model for the same 
orientation angles as in Figure 10. Independent of geomet- 
ric arrangement, spectra show a decrease in the polarization 
fraction with energy as they transition from unscattered to 
scattering-dominated. This transition occurs over the range 
2-4 keV, and involves a decrease in the net polarization frac- 
tion. 

Although this decrease is generic, the precise value of the 
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Fig. 18. — Phase-averaged intensity (top), polarization fraction (middle), 
and polarization angle (bottom) as a function of photon energy for models 
with different A(/>ims. Magnetic and viewing orientations con'esponding to 
drot = 45° and 9\as = 70° are assumed. Curves have the same meaning as in 
Figure 13. 



polarization fraction and angle are orientation dependent. A 
comparison of the middle panels of Figures 10 and 17 shows 
that the phase-averaged polarization fraction is significantly 
more dependent on orientation than the corresponding phase 
curves. This is largely due to the phase dependence of the 
polarization angle, and is most clearly illustrated by the solid 
black (6'rot = 45°, 9io, = 70°) and dotted red curves (9,ot = 70°, 
^los = 45°) in the bottom panels of Figure 10. In general, for 
^los < ^^rot, the polarization angle tends to sweep out the full 
180° range as the magnetic axis rotates on the plane of the 
sky. A larger polarization angle variation over a rotational pe- 
riod typically leads to a weaker phase averaged polarization 
fraction because the Stokes Q and U parameters have differ- 
ent signs at different rotational phases and tend to cancel. In- 
deed, the orthogonal rotator has the highest polarization frac- 
tion because the polarization angle shows the weakest phase 
dependence. 

The variation of the phase-averaged polarization angle with 
energy is also orientation dependent. In practice, it is difficult 
to extract useful constraints from this quantity alone, since it 
is the phase dependence of the polarization vector, both an- 
gle and magnitude, that determines the phase-average. Rather 
different orientations can provide very similar rotation of the 
polarization angle with energy, while rather similar geome- 
tries can give significantly different results. 

Figure 1 8 shows the variation of the phase-averaged observ- 
ables with twist angle for the same orientation as in Figure 13. 
Note that there is some diversity in the shape of the spectra for 
the total intensity, even though we vary 7,nax simultaneously 
with A0NS to keep the high energy spectral slope approxi- 
mately constant. As in the phase-resolved case, the polariza- 
tion fraction decreases with increasing A0ns, but Figure 18 
demonstrates that a drop in polarization fraction occurs for all 
photon energies. For increasing energy, the polarization an- 
gle rotates by an amount which increases monotonically with 
A(/)Ns- In principle, this dependence could be used to probe 
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Fig. 19. — Phase-averaged intensity (top), polarization fraction (middle), 
and polarization angle (bottom) as a function of photon energy for models 
with different seed photon polarization. Magnetic and viewing orientations 
corresponding to drot = 45° and 9\os = 70° are assumed. Curves have the 
same meaning as in Figure 15. 



the twist angle, but the rotation of the polarization angle de- 
pends on the orientation as well, so this needs to be indepen- 
dently constrained. 

Figure 19 compares the phase-averaged observables for dif- 
ferent seed photon distributions, assuming the same orienta- 
tion as in Figure 15. The curves correspond again to linear 
combinations of models tl0g20c90Eu and tl0g20c90Ou. 
At low energies (< 2 keV), the unpolarized and partially po- 
larized seed photon distributions have low polarizations frac- 
tions, as expected. The pure O-mode distribution has a lower 
linear polarization than the pure E-mode, because the small 
number of scattered photons at low energy with a net E-mode 
polarization which cancels some of the O-mode seeds. At 
higher energies (above ^ 5 keV), which are scattering domi- 
nated, mode exchange leads to a common value for the polar- 
ization fraction, independent of the seed distribution. 

The variation of the polarization in the transition region 
from 2-5 keV is very sensitive to the seed photon distribution. 
As the tl0g20c90Ou model transitions from being domi- 
nated by O-mode seed photons to E-mode scattered photons, 
there is a clear drop and subsequent rise in the polarization, 
which coincides with a 90° rotation of the polarization. A 
similar trend is seen in the partially polarized, but O-mode 
dominated (75% O-mode, 25% E-mode) model. In contrast, 
the E-mode dominated models show a continuous decline, 
with only modest rotation of the polarization angle. Hence, 
observations of the polarization fraction as a function of pho- 
ton energy in this transition region can discriminate clearly 
between predominately O-mode and predominantly E-mode 
seed photons distributions. This would help elucidate the 
mechanism providing the dominant contribution to the ther- 
mal luminosity. 

Figure 20 shows the effects of confining emission to po- 
lar caps. At high energies, the models all approach a com- 
mon value for the polarization fraction. The differences in 
the fraction of scattered photons primarily drive the observed 



17 





120 
115 
110 
105 
100 
95 
90 





E (keV) 

Fig. 20. — Phase-averaged intensity (top), polarization fraction (middle), 
and polarization angle (bottom) as a function of photon energy for various as- 
sumptions about NS surface emission regions and photon propagation. Mag- 
netic and viewing orientations corresponding to Smt =45° and 9[as = 70° are 
assumed. Curves have the same meaning as in Figure 16. 



differences in the polarization fraction at low energies. The 
highly polarized seed photons tend to escape more easily 
along the poles, particularly when the effects of light bend- 
ing are neglected. The exception is the model with emis- 
sion from only the south pole, because the observer primarily 
sees scattered photons for the assumed orientation (6'iot = 70°, 
^los = 45°). This also explains the large difference in the total 
intensity spectrum. An observer in the southern hemisphere 
(e.g. 6*108 = -45°) would see spectra that are qualitatively sim- 
ilar to the other polar cap models. In contrast to the phase 
curves of the polarization angle, which were identical for the 
different models, the phase-averaged polarization angles dif- 
fer due to the variation in the strength of polarization with 
rotational phase. 

5. IMPLICATIONS FOR GEMS 

We now perform simple estimates on the integration times 
required to obtain polarization information from magnetar 
candidates, based on current GEMS sensitivity parameters. 
The product of the GEMS effective area and modulation fac- 
tor rises from zero at an energy of 2 keV to its maximum 
around 3 keV, and then decreases steadily, in such a way that 
by 6 keV it is smaller by a factor of three from the peak 
[T. Kallman, private communication; see also Swank et al. 
(2010)]. This and the typical form of magnetar candidate 
spectra implies that the bulk of the measured polarization sig- 
nal will come from the transition between thermal and non- 
thermal parts of the spectrum. In §4 we showed that the po- 
larization fraction in this region can have large variations de- 
pending on the polarization level of the surface emission (and 
the temperature of the blackbody, which we do not change in 
this study). According to our results, the non-thermal tail is 
expected to produce a more uniform and predictable polariza- 
tion level, albeit at a significantly small flux level. Thus we 
base our estimates on trying to obtain a statistically significant 
measurement in this part of the spectrum, as it is the limiting 
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Following Weisskopf et al. (2010), we estimate uncertain- 
ties by assuming that the polarization noise follows a Pois- 
son distribution. The 1-sigma uncertainty in the modulation 
amplitude measured by the detector when rotated around the 
Une-of-sight is given by 
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(46) 



C=0.68 



where is the number of counts, A^4 = A^/IO"*, and C is the 
confidence level. The modulation amplitude produced by the 
source is 

ao = '«!!/,, (47) 

where m is the modulation of the detector, and XI^ is the linear 
polarization fraction (eq. [40]). We assume here that the con- 
tribution from the background is not important. Combining 
equations (46) and (47) allows determination of the number 
of counts needed to measure gq to within a given number of 
sigmas n^, 
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The 68% confidence contour of measured polarization angle 
corresponding to Aaia is (Weisskopf et al. 2010) 

ao sin xo + Aa i „ sin 



tanxiCT = ■ 



(49) 



flo cos Xo + Afl 1 o- cos ?/; 

where x xo ^re the measured and true polarization 
angle (eq. [41]), respectively, and ip is a free parameter. 
Finding the extrema of equation (49) yields tpexi = Xo ± 
arccos(-Aaicr/fl()). Setting xo = 0, to which the extrema of 
equation (49) are weakly sensitive, and keeping terms linear 
in Aflicr/flo yields 

AxifT = X-Xo = arctan ( — — ] — — . (50) 
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One can then straightforwardly obtain the number of counts 
required for a given uncertainty in the polarization angle mea- 
surement 

2 /r. 2 X 2 

(51) 
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In this simple analysis, the ratio N^/Nn oc {n„/ Ax)^ is insen- 
sitive to the modulation of the detector and the polarization 
fraction of the source. 

For a given magnetar candidate, the count rate in a given 
energy band is 

Ne^.e, = ^ ' (^^) Aeff(£)d£, (52) 

where [£[,£2] are the limits of the energy band, dN/dE the 
photon number flux, and Aeir is the effective area of GEMS. 
We obtain dN/dE from published fits to a handful of bright 
magnetar candidates in the sample of Kaspi & Boydstun 
(2010), which are selected to be in the quiescent stage'"^. It 

We reconstruct the number flux dN/dE using the published fit parame- 
ters in XSPEC for the best-fit blackbody plus power law models with neu- 
tral absorption (PHABS) included. In cases where model normalizations 
were not reported, we chose values to match the reported (energy band in- 
tegrated) fluxes or flux ratios. Since these models provide a reasonable fit for 
all sources under consideration here, they are sufficient for our approximate 
estimates. 
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TABLE 2 

Estimated GEMS Count Rates and Observational Times for Selected Magnetar 

Candidails 



Object iV2-io(s-') iV2^(s-i) iV4-io(s-i) ^nj^ (ks) Tl_ig(ks) Ref.' 



4U 0142+61 


3.27 


3.12 


0.15 


27 


101 


1 


IRXS J170849.0-400910 


0.767 


0.684 


0.083 


48 


180 


2 


IE 1841-045 


0.557 


0.475 


0.082 


49 


184 


3 


IE 2259+586 


0.723 


0.700 


0.023 


174 


653 


4 


SGR 1900+14 


0.113 


0.095 


0.018 


222 


833 


5 



"References: (1) White et al. (1996), (2) Rea et al. (2007b), (3) Moiii et al. (2003), (4) Woods et al. 
(2004), (5) Mereghetti et al. (2006) 



is well known that some of these sources experience long 
term changes in their flux and spectral properties (e.g., Kaspi 
2007), hence these estimates serve as a rough reference rather 
than a detailed prediction. 

The resulting count rates in photons per second are shown 
in Table 2. Alongside the number of counts over the entire 
sensitivity range of GEMS (2-10 keV), we show the rates 
over the soft and hard bands employed in §4.2 (2-4 keV and 
4—10 keV, respectively), which contain the transition between 
thermal and non-thermal components, and the pure power-law 
regime in most sources, respectively. We also show fiducial 
integration times for the 4 - 10 keV frequency band, 



'4-10" 



fX _ 
-'4-10" 



^4-lo 
W4-10' 



(53) 
(54) 



with m = 0.5, 11 = 20%, n„ = 3 and Axkt = 10 deg. Scal- 
ing the numbers with respect to this choice of parameters is 
straightforward. 

Obtaining just one phase-averaged data point in the pure 
power-law region of the spectrum requires integration times 
of the order of tens of kiloseconds or more. Because the ratio 
of count rates in the 2-4 keV to 4- 10 keV bands is 5-30, 
one would be guaranteed a higher energy resolution below 
4 keV, unless the seed photons are very unpolarized. Even 
then, one would be able to detect a single data point with a 
4.4% polarization fraction at 3-a in the band 2-4 keV for 
4U 0142+61. An unpolarized source of seed photons would 
require a cancellation of the contribution from returning cur- 
rents and deep cooling (§2.2), which we consider to be a 
less likely scenario. Thus prospects for obtaining indepen- 
dent clues about the mechanism producing the non-thermal 
X-rays are certainly feasible for several of the brightest AXPs 
and SGRs. 

The observational requirements for phase-resolved po- 
larimetry, on the other hand, increase linearly with the number 
of data points per rotation period, if the same accuracy wants 
to be maintained. Measurement of the polarization fraction 
with 5 points in phase is not very demanding, at least for the 
brightest objects in Table 2. The Icr uncertainty in the polar- 
ization angle using T^^q per phase point is Ax\a — 20 de- 
grees, which is satisfactory if the spatial scale of the twist is 
large and the orientation is favorable (e.g.. Figures 10, 13, and 
15). Accurately capturing smaller variations in polarization 
angle in both spatial scale and amplitude (e.g., lowest panel 
in Figure 12), would become much more demanding, with 
count rates easily exceeding 10^ s even for 4U 0142+-61. 



A measure of the polarization angle would be especially 
useful around outbursts, as it has already been seen that 
sources experience persistent changes in torque and pulse pro- 
files (e.g.. Woods et al. 2004, 2007). Simultaneous monitor- 
ing with RXTE and high resolution spectra with Chandra or 
XMM-Newton would be the ideal complement to deep polar- 
ization observations. 

6. SUMMARY AND DISCUSSION 

We have studied the effects of resonant cyclotron scattering 
in a twisted magnetosphere on the outgoing polarization 
signal of magnetars. We have employed simple prescriptions 
for the magnetic field geometry, particle energy distribution, 
and seed spectrum. When the pair multiplicity is not too 
large, the dielectric properties of the magnetosphere are 
dominated by vacuum polarization out to large distances 
(eq. [17]), with the following implications: 

1. - There is an effective separation of the regions where 
resonant cyclotron scattering and polarization eigermiode 
coupling take place, with the consequence that scattering 
occurs mostly in the adiabatic limit (Fig. 2). Because the 
light cylinder is safely outside of the zone where polarization 
modes freeze, circular polarization is always less than 10% 
(e.g.. Figure 4). 

2. - Regardless of the polarization distribution of surface 
photons, resonant cyclotron scattering by transrelativistic 
charges introduces a characteristic polarization fraction 
of the order of ~ 10-30%. This number arises from the 
mode exchange probabiUty, dependence of optical depth on 
polarization, multiple scattering, and geometry of the system 
(§4.1). The transition from surface-dominated polarization to 
scattering-dominated takes place over a region of the spec- 
trum extending over the interval ~ 2 — 5 keV for a thermal 
spectrum with kT = 0.4 keV (Figures 17-20). Measurement of 
phase-averaged polarization fraction allows an unambiguous 
identification of resonant comptonization as the mechanism 
producing the non-thermal soft X-ray emission. This is 
particularly relevant for sources such as XTE J1810-197, 
whose spectrum can be fit with a comptonized model or 
two-temperature blackbody (Halpem & Gotthelf 2005). 

3. Phase averaged measurements can also discern the domi- 
nant polarization mode of seed photons (Fig. 19), allowing 
an independent constraint on the mechanism providing the 
thermal luminosity. A dominance of E-mode photons in the 
thermal spectrum would favor deep cooling of the neutron 
star (Thompson & Duncan 1996), setting constraints on 
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magnetized atmosphere models (e.g., van Adelsberg & Lai 
2006). An O-mode dominance in the unscattered component 
would indicate thermal emission from returning currents 
hitting the stellar surface (Thompson et al. 2002). 

4. - Measurement of the phase-resolved polarization angle 
allows measurement of the spatial extent of the twist angle 
at radii ^ 1007?ns- Because the mode exchange probability 
always favors the E-mode regardless of system parameters 
(eq. [16]), scattered photons are predominantly in this linear 
eigenstate and their adiabatic propagation keeps track of 
the magnetic field geometry, with little sensitivity to input 
parameters. Photons remain mostly linearly polarized as 
eigenmodes couple, a direct consequence of the dominance of 
vacuum polarization and the slow rotation of magnetars. For 
each photon, the polarization vector preserves its orientation 
relative to the direction of the magnetic field at the point of 
freezing. In the case of a globally twisted dipole, it may even 
be possible to obtain the geometric parameters of the system 
by inverting the polarization angle as a function of phase. 
We have not explored this possibihty (e.g., lower panels of 
Figures 10 and 14). It is worth keeping in mind that this 
effect is dependent on magnetospheric scattering only in that 
a particular eigenmode is preferred: in principle it is also 
present in the unscattered component (e.g., Heyl & Shaviv 
2000). 

5. - The asymmetries in the angular distribution of scattered 
radiation are due to the combination of charge motion along 
field Unes and relativistic aberration. This also imparts 
an asymmetry in the polarization fraction. If scattering 
charges consist predominantly of electrons moving from one 
magnetic pole to the other, then a north-south asymmetry 
is introduced in the polarization fraction. Phase-resolved 
polarization fraction wiU display a peak displaced by 
one-half of a phase relative to the peak in intensity if the 
particle energy distribution is uni-directional (Figure 12). 
If the magnetosphere is dominated by pairs at these radii 
(^ IO/Jns)^ then two symmetric peaks are expected in both 
the polarization distribution and the radiation flux, barring a 
strong asymmetry in the surface emission. 

6. - Based on the expected sensitivity of GEMS and the 
spectra of a few bright magnetars, we estimate that phase- 
averaged polarimetry requires several tens to a few hundred 
kilosecond exposures per source for obtaining at least one 
data point in the non-thermal tail of the spectrum (Table 2). 
Phase-resolved polarimetry increases the requirements pro- 
portional to the number of points in phase for fixed accuracy. 
Detailed measurements of changes in the polarization angle 
as a function of phase seem feasible only for the brightest 
sources, and would require integration times approaching or 
exceeding a 10^ s. 

Our analysis ignores the contribution of the magnetospheric 
plasma to the dielectric tensor. This is a good approxima- 
tion if the density of charge carriers satisfies He ~ J/{ec), as 
found in Beloborodov & Thompson (2007). That study, how- 
ever, focused on the inner part of the closed magnetosphere 
at radii r ~ /Jns> where a number of effects present at larger 
radii are absent. Pair cascades in the closed field circuit could 
be triggered by energy deposition from current driven insta- 
bilities close to the star (Thompson 2008a). A significant 
number of pairs at IQQRjis then alters the relative contribu- 



tion of plasma and vacuum to the dielectric tensor (eq. [17]). 
When the plasma component is non-negligible, and a; ^ cjc, 
propagation eigenmodes become increasingly circular as the 
charge density is increased, for most propagation directions 
(e.g., Meszaros 1992). Capturing this effect, particularly if 
the transition from vacuum to plasma dominated modes oc- 
curs around rpi, requires a more careful consideration of the 
distribution of currents in the magnetosphere than our sim- 
ple space-independent particle energy distribution. We have 
opted for a self-consistent (but potentiaUy incomplete) first 
step in modeling this problem. 

We have also ignored the contribution of scattering by ions 
close to the star (e.g., Thompson et al. 2002). Since the di- 
electric properties of the plasma would not change relative to 
the case analyzed here, we expect only quantitative changes 
relative to scattering by uni-directional electrons. Ion scat- 
tering would still be in the adiabatic limit, and modes would 
freeze at the same distance from the star as in the case of pho- 
tons scattered by electrons. Hence, we would not expect sig- 
nificant changes in the phase dependence of the polarization 
angle. Changes would mostly affect the intensity and polar- 
ization fraction, and would be due to general relativistic ef- 
fects (in photon propagation and magnetic field, the second 
of which we are ignoring) as well as due to reflection or ab- 
sorption by the star, whose radius would now be comparable 
to the resonant surface. And obviously photons that escape 
the ion resonance surface can still be subject to scattering by 
electrons at larger radii, so the net observed anisotropics could 
increase. 

Most of our calculations assume the simplest possible pre- 
scription for the seed photons: a blackbody of uniform sur- 
face temperature, 100% polarized, and isotropically emitting, 
in order to isolate the magnetospheric imprint on the polariza- 
tion signal. Out of these assumptions, we expect those con- 
cerning the surface temperature and beaming distribution to 
potentially alter our results once they are relaxed. To gauge 
the effects of a non-uniform surface temperature, we have 
explored the extreme case of narrow polar caps around one 
or both magnetic poles, motivated by results from thermal 
conduction calculations with a strong internal toroidal field 
(Perez-Azorin et al. 2006). We find minor differences relative 
to the uniform surface temperature distribution unless one po- 
lar cap emits substantially more than the other (Figs. [16] and 
[20]). A radiation intensity that has a non-trivial distribution 
around the local magnetic field (e.g., van Adelsberg & Lai 
2006) could also distort the signature due to the directionality 
of particles. We have not addressed this issue here, as it would 
require a realistic atmosphere for self-consistency. 

Given our assumptions, the polarization angle is an almost 
perfect tracer of the magnetic field geometry at ~ 100/?ns 
from the star. In this study we have employed a very sim- 
ple field geometry, which is self-similar and thus having a 
local direction independent of radial distance from the star. 
In a more realistic magnetosphere, the twist angle is not nec- 
essarily global, and evolves both in time and space (Thomp- 
son et al. 2002; Beloborodov & Thompson 2007; Thompson 
2008a; Beloborodov 2009). Measurements of the magneto- 
spheric twist at large radii would be most useful during and 
following periods of burst activity, where the magnetosphere 
is expected to undergo a significant rearrangement, as shown 
by changes in the pulse profiles before and after outbursts 
(e.g., Woods et al. 2004, 2007). A pulse profile with strong 
asymmetries and small-scale structure, if associated with lo- 
calized twists, should be concurrent with a phase-resolved 
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Fig. a 1 . — Coordinate systems and geometric transformations employed in the Monte Carlo calculation, labeled by roman numbers. The vectors M, Q, and O 
denote the magnetic axis, rotation axis, and line-of-sight (see also eqns. [42] and [43]). The arable numerals indicate the order in which each rotation matrix is 
applied (see text). 



polarization angle that displays high temporal frequency and 
small amplitude components. 
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APPENDIX 

A. COORDINATE TRANSFORMATIONS 

We present here all the coordinate transformations employed in the calculation. We denote the coordinate system with the 
magnetic axis M in the z-direction by the subscript M. Each alternative coordinate system employed is labeled by roman number 
subscripts on vectors. Each of them is illustrated in Figure Al and explained below. For conciseness, we write the transformations 
in terms of the elementary rotation matrices 



1 

cost? 
smt9 





-sint? 

C0Sl9 



COS I? sin I? 

1 
-sini? cos6 



cos?? 
sin?? 




-sin?? 
cos?? 
1 



(Al) 



which when appUed to a column vector x produce a counterclockwise rotation by an angle ?? around each of the cartesian 
coordinate axes, respectively. We also make use of the more general angle-axis formula: given a unit vector u, the rotation matrix 
that produces a counterclockwise rotation by an angle •& around u = UxX+ Uyy + u,z is 

Mj + (1 -M^)cos?? UxUy{l-cos^)-UzSm^ UxU^il- cost)) + Uy sin 

/?„(??)= M^Mj,(l -COS??) + MjSini? My + (l-My)C0S'!? Mj,Mj(l - COS??) - M^sin?? 

UxU^{l-cos'&)-UySin'& UyU^{l-cos'&) + UxSm'& m| + (1-m|)cos?? 
with (Al) corresponding tou = x, y, and z, respectively. 



(A2) 



A.l. Light Bending (I) 

Since a photon orbit in a Schwarzschild metric has two conserved quantities (energy and angular momentum), one can map 
it to the ^ = tt/2 plane in spherical polar coordinates without loss of generality (e.g., Schutz 2009). We choose the x-axis to be 
coincident with the position vector f, and the z axis along fx k. The photon wave vector is transformed according to 

ki = R,{-ip)R,{'&)kM. (A3) 



21 



The first rotation around u = (f xk) x M /\f x k\, with cos'd = {fxk)-M/\fxk\, puts the traj ectory on to the equatorial plane. We 
then rotate the radial vector f towards x, with 

The normalized wave vector evolved in equations (29)-(32) corresponds to the components of kj in spherical polar coordinates. 
The transformation is ill-defined for radial propagation, which is nonetheless unaffected by light bending. In practice, we revert 
to planar propagation whenever 1 - |A: • r| < 10"^. 

A.2. Polarization Evolution (II) 

We start each integration of equations (21)-(22) in a frame such that the wave vector k points in the z direction and the magnetic 
field lies in the x-z plane. The transformation for the electric field vector A is 

A// = Rz(-<pB)Ry{-ek)Rz(.-<Pk) Am, (A5) 
where and are the polar and azimuthal angle of the wave vector relative to the magnetic axis, and 

[RA-OijR,(-ok)B]y 

tan (pB = — '■ ^ . (A6) 

[Ryi-9,)R,i-A)Bh 

A3. Rotational Phase (III) 

Our coordinate system for phase-dependent quantities is such that the rotational axis points in the z-direction, and at phase 
zero both the magnetic axis M and the direction to the line-of-sight O are in the x-z plane, at angles ^rot and ^los from Cl, 
respectively. What is needed are the coordinates of the hne-of-sight O relative to the magnetic axis, so that stored quantities are 
retrieved to construct phase curves. The transformation is 

Om = Ryi-erodRz(-<Pm) dm, (A7) 

with Oni = smOiosXiii + cosOiosZiii, and c^^rot the rotational phase. 

A.4. Stokes Parameters (IV) 

A rotation of the polarization plane around the wave vector affects only Q and U (Chandrasekhar 1960). Counterclockwise 
rotation by an angle (p yields 

Q' = Qcos2ip + Usm2ip (A8) 
U' = -Qsm2ip + Ucos2ip. (A9) 

For storage after photon escape to infinity, we take the x axis in the ^-M plane. For phase dependent quantities, however we 
rotate this axis to coincide with the k-Cl plane. The needed angle ip is obtained by rotating Om so that it lies along the z axis, 
and then computing the azimuthal angle between -Mjv and -Cljv, thus 

tmin -ip)= - y - , (AlO) 

where Mm =Rz((t>wt)RyiOm)z (eq. [A7]). 
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